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In the wake of previous studies on the rattling-and-jumping diffusion in smectic liquid crystal 
phases of colloidal rods, we analyze here for the first time the heterogeneous dynamics in columnar 
phases. More specifically, we perform computer simulations to investigate the relaxation dynamics 
of a binary mixture of perfectly aligned hard spherocylinders. We detect that the columnar arrange- 
ment of the system produces free-energy barriers the particles should overcome to jump from one 
column to another, thus determining a hopping-type diffusion. This phenomenon accounts for the 
non-Gaussian inter-column diffusion and shows a two-step structural relaxation which is remarkably 
analogous to that of out-of-equilibrium glass-forming systems and gels. Surprisingly enough, slight 
deviations from the behavior of simple liquids due to transient cages is also observed in the direction 
perpendicular to this plane, where the system is usually referred to as liquid-like. 



I. INTRODUCTION 

Liquid crystals (LCs) are phases of matter in which 
the anisotropy of the particles determines under specific 
conditions a partial spontaneous breaking of the spa- 
tial symmetries of the system, thus manifesting features 
in between the crystalline solid and the isotropic liquid 
phase. The notion that entropic effects alone are suffi- 
cient to drive the self-assembly of ordered liquid crystal 
phases is well established in colloid science drill- Model 
systems of hard particles constitute the natural choice 
to describe the phase and aggregation behavior of most 
colloidal systems, as the main interactions established 
between their particles have a repulsive, steric origin. As 
a consequence, a hard-particle fluid does not have inter- 
nal energy and minimizing its free energy is equivalent 
to maximizing its entropy. 

In his seminal work, Onsager showed that mere hard- 
core repulsions between infinitely thin rigid cylinders 
are able to determine an entropy-driven phase transition 
from the isotropic to the nematic phase, and hence the 
existence of a spontaneous oricntational order [l[. The 
evolution of simulation techniques in the past decades 
allowed to investigate the more realistic case of rods 
with finite size, showing that by varying the length-to- 
diameter ratio one- (smectic) , two- (columnar) and three- 
crystal) dimensional translational ordered phases can 
be encountered [f| 0|- Further more accurate studies 
showed that for a monodisperse system of both aligned 
[H and freely rotating Q hard spherocylinders the colum- 
nar phase happens to be metastable with respect to the 
smectic phase for each value of the length-to-diamcter 
ratio. The complexity of the phase behavior of linear 
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particles becomes even more pronounced by proceeding 
from monodisperse systems to mixtures [10l - ll4l |. In par- 
ticular, size polydispcrsity introduces a sensible change 
in the phase behavior of rod-like particles as their pack- 
ing is not as effective as that of monodisperse rods. In 
fact, it was observed that in a system of hard sphero- 
cylinders the formation of smectic layers can be inhibited 
by introducing a length bidispersity, in such a way that 
the columnar phase can become thcrmodynamically sta- 
ble [T5L ITlj . Entropy-driven columnar phase transitions 
have been observed in monodisperse systems of disc-like 
particles, such as cut spheres or oblate spherocylinders 
[15l - [l8f . As far as rod- like particles are concerned, theo- 
retical studies indicated that the columnar order can be 
observed not only in bidisperse mixtures, but also in more 
realistic polydisperse systems of parallel cylinders fl9| . 
Polydispcrsity is not the only element which favors the 
stabilization of the columnar phase in a system of rods, 
as was shown in Ref. [2^] where a monodisperse system 
of soft-core rods was considered. On the other hand, the 
effect of rod flexibility in stabilizing the columnar phase 
is still under debate [2l| - |24| . 

With the improvement in understanding and describ- 
ing static and thermodynamic properties of LCs, the in- 
terest towards the dynamics correspondingly increased. 
In particular, most of the studies in this direction were 
devoted to analyze the anisotropy of the diffusion by the 
measurement of the self-diffusion coefficients in differ- 
ent mcsophases p5l - l28j . For lyotropic LCs these results 
found good overall agreement with experiments based 
on techniques such as fluorescence recovery after photo- 
bleaching (FRAP) Hi-ill. On the other hand, only few 
studies focused on the analysis of the dynamical phe- 
nomena at the single-particle level, where fluorescence 
microscopy was applied to investigate the LC phases in 
colloidal suspensions of fd virus [3J, H3] • This approach 
allowed to observe for the first time the mechanism of 
intcrlayer diffusion, or permeation, which characterizes 
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the dynamics of the system in the smectic phase [34J. 
In fact, the layered structure of smectic LCs, which de- 
termines an effective periodic mean-field potential, influ- 
ences the motion in the direction perpendicular to the 
smectic layers with the appearance of jumps of the order 
of the rod length. As a result of this quasi-quantized dif- 
fusion, at intermediate time intervals one can distinguish 
between "slow" particles, which rattle around the center 
of a column, and "fast" particles, which jump to another 
column. This heterogeneous dynamics accounts for devi- 
ations from Gaussianity in the diffusion, a phenomenon 
which is also well known in homogeneous complex fluids 
such as supercooled liquids [H, ® an d gels H3 ■ In 
the light of these results, a theoretical approach based on 
dynamic density functional theory using the second virial 
approximation focused on the role played by the local 
fluid structure of the system, which competes with the 
permanent barriers due to the smectic structure, thus in- 
creasing even more the analo gies with fluids close to a dy- 
namical arrest transitio n I39l.l40j ] . Simulations on aligned 
[4lj and freely rotating (42T 43 1 hard spherocylinders con- 
firmed qualitatively both these studies and pointed out 
the effect of dynamical heterogeneities on the structural 
relaxation of the system, which deviates from the expo- 
nential decay exp ected for simple fluids. Furthermore, 
in Refs. 
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dynamics in the smectic phase of hard spherocylinders is 
tightly related to cooperative motion of particles through 
the smectic layers. 

Following this line of research and motivated by a re- 
cent experiment on fd virus [44[, we investigate in this 
paper the dynamics of a binary mixture of rod-like parti- 
cles which exhibits a stable columnar liquid crystal phase. 
Using Monte Carlo (MC) simulations we are able to study 
for the first time the dynamical heterogeneities which 
arise from the columnar structure and their effect on the 
diffusion and on the long-time structural relaxation of 
the system. Furthermore, we measure the height of the 
inter-column energetic barriers and compare our results 
with the simulation data available for the smectic phase 
[4lTj43j . Conventional MC dynamics has been extensively 
applied to study the ori gin of non-exponential relaxation 
of glass-forming liquids [45l - l49j . MC and Molecular Dy- 
namics (MD) simulations of Lennard- Jones fluids [501 ] 
showed that the particle dynamics at large time scales 
are the same, but are different at small time scales, as 
the stochastic motion of particles cannot be detected by 
the deterministic approach of MD simulations. 

The paper is organized as follows. In section II we in- 
troduce the model, the simulation details and the physi- 
cal properties which were measured in order to describe 
the dynamics of the system. The results of these mea- 
surements, which confirm the presence of dynamical het- 
erogeneities also in the columnar phase and are in general 
agreement with the observations in the smectic phase, are 
discussed in section III, whereas in section IV we illus- 
trate our conclusions. 



II. MODEL AND SIMULATIONS 

We study a system containing N = 1600 perfectly 
aligned hard spherocylinders with aspect ratio L* = 
L/ D, where L and D arc, respectively, the length and di- 
ameter of a cylindrical body capped by two hemispheres 
with diameter D. The phase diagram of a monodisperse 
system containing such rod-like particles shows stable ne- 
matic, smectic, and crystal phases, but lacks a stable 
columnar phase in the range < L* < oo Q. Stroobants 
studied the phase behavior of bidisperse systems of hard 
rods and found that the bidispcrsity can favor and sta- 
bilize columnar order over smectic order [Io| . Therefore, 
to prevent the formation of smectic layers, we investigate 
a binary mixture of hard spherocylinders with the same 
diameter D (used as our unit of length), but different 
lengths L\ and L 2 , with L\ > L 2 . In this model, the 
rotational degrees of freedom are frozen out and hence 
the particles are forced to be aligned along a common 
nematic director, oriented along the z axis. The relative 
concentration of the two species is set in such a way that 
the binary mixture is kept at its equivalence point, where 
the volume fractions of each component are the same. 
The phase diagram at fixed L 2 = 1.0 displays a region 
of stability of the columnar phase which increases with 
L\ and disappears at L?< 1.6, where a nematic- smectic 
transition is observed [101 ] . Here we study a columnar or- 
dered binary mixture of rods with L\ = 2.1 and L 2 = 1.0, 
and relative concentrations x\ = N\/N = 0.375 and 
x 2 = N 2 /N = 0.625, respectively. For lower pressures, 
this columnar phase transforms into a nematic phase, 
while for higher pressures it freezes into a crystal phase. 

We performed standard MC simulations in a rectangu- 
lar box of volume V with periodic boundary conditions. 
To equilibrate the columnar phase, we performed runs in 
the isobaric-isothermal (NPT) ensemble, where the par- 
ticle moves were accepted according to the Metropolis al- 
gorithm [5l|, that is if no particle overlap was detected. 
Each MC cycle consisted of N attempts to displace a ran- 
domly selected particle, plus an attempt to modify the 
box volume with independent changes of the three box 
sides. The system was considered to be in equilibrium 
when the volume reached a stationary value within the 
statistical fluctuations. We run simulations at several re- 
duced pressures P* = /3PD 3 , where fi = l/fc^T, ks is 
the Boltzmann's constant and T the absolute tempera- 
ture. In particular, we equilibrated a nematic phase at 
P* = 2.5 (packing fraction i] = N(x\V\ + x 2 V2)/V = 
0.470 with Vi (i = 1, 2) the single particle volume) which 
is very close to the nematic-columnar transition, and 
three different columnar phases at P* = 3.0 (rj — 0.535), 
3.5 (rj = 0.563) and 4.0 (77 = 0.580). In all these cases, our 
starting configuration consisted of a highly packed colum- 
nar structure with the rods randomly located along the z 
direction, and hexagonally ordered in the xy plane. The 
minimum number of MC cycles needed for an equilibra- 
tion run was 5 x 10 5 , and was followed by a production 
run of 2 x 10 6 MC cycles in the canonical (NVT) en- 
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scmblc to simulate the relaxation dynamics and evaluate 
all the physical properties of interest. In this case, the 
box volume was kept fixed to prevent unphysical collec- 
tive moves which do not mimic the Brownian dynamics 
of the particles properly. It was proved that in rod sus- 
pensions the contribution of hydrodynamics can at first 
approximation be neglected with respect to steric effects, 
which result in excluded volume interactions [52j . Under 
these conditions, the MC approach offers an important 
tool to study the dynamics of colloids, since, in spite of 
its intrinsically non-dynamical nature, it is able to re- 
produce the Brownian diffusion typical of such kind of 
systems (53|. To pursue this goal, one must set a small 
enough maximum MC displacement, typically of the or- 
der of one tenth of the shortest dimension of the particle. 
The optimal value of the mean particle displacement is 
strictly linked to the acceptance rate and hence to the 
CPU time per simulation run. If the displacement is 
chosen too small, the system would need longer runs to 
properly explore the configurational space. The same ef- 
fect is expected with too big displacements, as most of 
the particle moves would cause overlaps and therefore 
would be rejected. However, any reasonable and con- 
venient choice of the average step size should not affect 
the dynamics at long time scales. We tested that, apart 
from an overall scaling, this was the case. The maxi- 
mum displacement was fixed to give an acceptance rate 
of roughly 50% per move. Furthermore, to take into ac- 
count the non-spherical shape of the particles, the max- 
imum MC displacement was set in such a way that at 
short times it reproduces the anisotropic diffusion of a 
single colloidal rod. More specifically, this means that 
the ratio between the maximum MC displacement in the 
xy plane Ax max = Ay max and in the z direction Az max 
must be set as 



Ax Ti 



A;. 



D 



(1) 



where we denoted with D± and D\\ the short-time self- 
diffusion coefficients of the rod in the direction paral- 
lel and perpendicular to its long axis, respectively. To 
give an estimate for the ratio D±/Dn, which solely de- 
pends on the geometry of the particle, we referred to the 
semi- empirical expression derived in Ref. (54|. In that 
work, the authors used a numerical approach to evaluate 
the translational self-diffusion coefficients of a cylindri- 
cal particle for different values of the length-to-diamctcr 
ratio p, finding good overall agreement with experimen- 
tal data in the range 2 < p < 30 [55[. A least-square 
quadratic fitting in p^ 1 of the data allowed the authors 
to give an expression for the transverse and longitudinal 
self-diffusion coefficients as functions of the parameter p. 
Since here we consider spherocylinders, we evaluated the 
ratio Dj_/Du by setting p = (L + D)/D, arguing that 
possible deviations given by the non-cylindrical shape of 
the particles were irrelevant. According to the above- 
mentioned expression, the ratio between the maximum 



MC displacement perpendicular and parallel to the z axis 
was set to 0.92 for particles of species 1 and 0.94 for those 
of species 2. Moreover, since the transverse section of 
the particles is the same for the two species, we set the 
same maximum MC displacement along the z axis for the 
two components. Once the short-time self-diffusion coef- 
ficients are known, it is possible to introduce a time scale 
defined by r = D 2 /D tr , where the total translational dif- 
fusion coefficient D tr = ((D\\) + 2(D±))/3 is evaluated in 
terms of the longitudinal and transverse short-time dif- 
fusion coefficients averaged over the two species . 

In order to analyze the heterogeneous diffusion and the 
structural relaxation of the system, the following physi- 
cal properties were calculated: (i) the transverse mean- 
field potential, (ii) the self-part of the van Hove function 
(SVHF), (iii) the distinct-part of the van Hove function 
(DVHF), (iv) the mean square displacement (MSD), (v) 
the non-Gaussian parameter (NGP) and (vi) the self-part 
of the intermediate scattering function (SISF). 

Transverse mean-field potential. In a liquid crystal 
phase characterized by columnar order the translational 
invariance is spontaneously broken in the plane perpen- 
dicular to the nematic director. This gives rise to a non 
homogeneous (relative) probability Tti(x,y) of finding a 
particle of species i = 1, 2 at position (x, y) in the plane 
perpendicular to the nematic director. The effective en- 
ergetic barrier which tends to confine the particle inside 
a column i s g iven by the mean- field potential Ui(x,y), 
defined as [34| 
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where the proportionality constant is chosen in such a 
way that the minima of the potential is set to zero. 

Self-part of the van Hove function. The heterogeneous 
dynamics and hopping-type inter-column diffusion can 
be quantitatively described by the SVHF [56^ 



which measures the probability distribution for a particle 
displacement r in a time interval t. Since the present sys- 
tem is characterized by a translational symmetry along 
the nematic director, it is natural to separately study 
the diffusion along the z axis and in the xy plane. This 
can be done by partially integrating the SVHF on the xy 
plane to get its longitudinal component 
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and along the z axis to get its transverse component, 
which can be further averaged over the azimuthal angle 
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In the above equations (Rj(t), Zj{i)) is the position 
of particle j at time t, 5 is the Dirac delta, (...) stands 
for an ensemble average and the index 2ir indicates an 
additional average over the polar angle, which defines 
the bidimensional vector R = (x, y) with modulus R = 
R|. It should be noticed that for freely diffusive particles 
these functions are described by a Gaussian. 

Distinct-part of the van Hove function. A description 
of the influence of the surrounding particles background 
on the single particle diffusion is given by the DVHF, 
which is the probability distribution on the relative po- 
sition r of two different particles at different times 



G d (r,t) = ±/ J2 5(r-r,-(t + to)+r i (t )) 
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In order to separately study the caging regime related to 
the local fluid structure in the direction longitudinal and 
transverse to the nematic director, we actually measured 
the DVHF partially integrated over the size occupied by 
the spherocylinder in the xy plane and along the z axis, 
defined, respectively, by 
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where we set L = min{I/i, L 2 } and 6 is the azimuthal 
angle of r in the xy plane. 

Non- Gaussian parameter. The deviations of the dif- 
fusion from Gaussian behavior can be estimated by the 
NGP, defined as [E3 



a 2 (t) = 
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where Ar(i) is the displacement of a particle during a 
time interval t. The parameter d corresponds to the num- 
ber of dimensions considered, so that d = I for the linear 
diffusion longitudinal to the nematic director (a 2tZ (t)) 
and d = 2 for the planar transverse diffusion {ot 2tXy (t)). 
As pointed out in Rcf. [46j], when treating mixtures one 
should be careful in not taking into account trivial non- 
Gaussianity due to a size-dependent particle mobility. In 
order to describe the dynamical heterogeneities related 
exclusively to the permanent barriers of the LC phase, 
one has to calculate the NGP a 2 as defined in Eq. ([9]) 



for each species i = 1, 2 separately and then perform an 
average weighted over the concentrations Xi, i.e. 



(02(f)) = x l a 2 1 \t)+x 2 a 2 2) {t). 
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With these definitions heterogeneous diffusion can be de- 
tected when the NGP deviates from zero value. 

Self-part of the intermediate scattering function. The 
structural relaxation of the system is conveniently de- 
scribed by measuring the self-part of the intermediate 
scattering function 
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which describes the density auto-correlations decay in 
the reciprocal space. Since most of the relevant struc- 
tural information is contained at the first peak k* of 
the structure factor, we can focus on the transverse and 
longitudinal relaxations by evaluating this function at 
(k*,k*,Q) and (0,0, fc*) respectively, so that F sxy (t) = 
F s ((k*,k* y ,0),t) and F s>z {t) = F s ((0, 0, k* z ),t). 



III. RESULTS 

When the difference in length between the two compo- 
nents of a binary system of aligned hard spherocylinders 
is sufficiently high, the system undergoes a transition 
from a nematic to a columnar phase by increasing the 
pressure. The structure of the columnar phase is charac- 
terized by the development of an hexagonal order in the 
plane perpendicular to the nematic director [To| . This is 
the case for our parameter choice L\ — 2.1 and L\ = 1.0 
as illustrated in Fig. [T] where two typical configurations 
in the nematic (P* — 2.5) and in the columnar phase 
(P* = 4.0) are compared. 

The development of a long-range translational order 
allows to interpret the diffusion as the motion of a single 
particle subject to a periodic mean-field potential U(x, y) 
as defined in Eq. (|2|). This approach was applied in 
experiments Hi| and simulations [4ll - [43j to character- 
ize the hopping-type diffusion in smectic liquid crystals 
along the nematic director. These authors found that the 
free-energy cost for the layer-to-layer diffusion is in the 
order of few ksT per particle, depending mostly on the 
packing of the system, but also on the anisotropy and ro- 
tational degrees of freedom of the rods. Since the system 
studied here is composed of two species, the mean-field 
potential was evaluated separately for each component 
and is shown in Fig. [5] for several pressures. The min- 
ima of the potential correspond to the lattice position, 
and the height of the energetic barriers gives a quantita- 
tive description of the energetic demand associated to a 
column-to-column jump. In order to estimate the height 
of the energetic barriers for each species as a function of 
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FIG. 1. (Color online) Side and top views of two typical 
configurations in the nematic (P* = 2.5, left figure) and in 
the columnar phase (P* — 4.0, right figure) of a binary mix- 
ture of perfectly aligned hard spherocylinders with length-to- 
diameter ratios L* — 2.1 and L% = 1.0 and relative concen- 
trations xi = 0.375 and X2 = 0.625. 




FIG. 2. (Color online) Mean-field effective potential U (x, y) in 
units of fcsT in the bulk columnar phase of a binary mixture of 
perfectly aligned hard spherocylinders at P* = 3.0, P* = 3.5 
and P* = 4.0 (from top to bottom). The images on the 
left correspond to the long rods (species 1), whereas those on 
the right to the short ones (species 2). In order to ease the 
visualization, the black lines at the top of each graph identify 
the isopotential points in the xy plane with increments of 
3fc s T. 



pressure, we report in Fig. [3] a transverse section of the 
energ y la ndscapes in Fig. [21 Following the procedure in 
Ref. [41(, the experimental points in Fig. [3] were fitted 
with a function 



U(R) = Y. Uk 
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with Uk and h fit parameters and n = 5. As expected, the 
height of the potential barrier increases with the packing 
fraction and with the particle anisotropy, as already ob- 
served in Ref. [H}. At significant packing fractions, we 
detect that the energetic barriers appear higher for long 
rods. In particular, at P* = 3.5 and 4.0 the column-to- 
column jumps become so rare that the associated statis- 
tics is too poor to furnish a precise estimate of the barrier 
height. In other words, the long rods are constrained to 
rattle in their cage, the jump to a neighboring column be- 
ing too demanding. This is due to the fact that at high 
packing fraction no MC configuration showed a long rod 
in the region between the columns, with the result that 



the mean-field potential was characterized by an unphys- 
ical divergence. Furthermore, one should notice that the 
typical height of the barriers, close to and even higher 
than lOfcffT 1 , is significantly higher than in the smectic 
phase |4l] - l43| . This can be seen by comparing our data 
at P* = 3.5 (rj = 0.563) with those in Ref. gH for the 
smectic phase of a system of aligned hard spherocylin- 
ders with L* = 5.0 at pressure P* = 5.0 (77 = 0.563). 
In the latter the height of the energy barrier reaches a 
value close to %ksT, which is expected to be even lower 
for shorter rods, as noticed in Ref. [42j. 

The effect of the periodic mean-field potential can be 
further appreciated in Fig. 01 where we show a typi- 
cal trajectory projected on the xy plane of a long and 
a short particle at P* = 3.0. The difference with the 
Gaussian diffusion typical of a simple liquid, where the 
particle trajectories resemble the behavior of a random 
walker, is evident. In this case, the dynamics of the sys- 
tem is characterized by a hopping-typc diffusion, in which 
each particle tends to rattle around the center of a col- 
umn until it finds suitable conditions to overcome the 
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FIG. 3. Transverse section of the mean-field effective potential 
in Fig. [5] for a binary mixture of long (a) and short (b) hard 
spherocylinders at P* = 3.0 (o), P* = 3.5 (A) and P* = 4.0 
(A). The solid lines are fits (see text). 
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FIG. 4. (Color online) A typical trajectory projected on the 
xy plane of a long (red line) and short (green line) rod in 
the columnar phase at P* — 3.0 after an interval of time 
At = 380t. 



energetic barrier and jump to another column in quasi- 
quantized steps. The spread in the total displacement 
between the two species is due to the higher barriers felt 
by the long rods, which significantly inhibit the inter- 
column diffusion. This behavior is observed in the whole 
range of pressures considered here, and its effects on the 
long-time relaxation dynamics of the system are crucial. 
More specifically, the long rods are expected to sample 




1 2 3 4 5 6 0.5 1 1.5 2 2.5 3 
R/D 



FIG. 5. Transverse component of the self-part of the van 
Hove function G 3 (R,t) as a function of R/D for a binary 
mixture of perfectly aligned hard spherocylinders with length- 
to-diameter ratio L* — 2.1 and = 1.0 and relative concen- 
trations xi = 0.375 and xi = 0.625 at t/r = 1 (solid lines), 
t/r = 10 (dashed lines) and t/r = 100 (dotted lines) for the 
system at pressure P* — 2.5 (a), P* = 3.0 (b), P* = 3.5 (c) 
and P* = 4.0 (d). 



the configurational space on a time scale which might be 
significantly longer than that needed for the small ones. 
As a consequence, the decay of the correlation functions 
is strongly affected by the slow diffusion of the long par- 
ticles, as we show later on. 

In order to quantitatively study the distribution of dis- 
placements of the particles at different pressures and after 
different time intervals, we report in Fig. [5]and[6]the self- 
part of the van Hove function (SVHF) in its transverse 
and longitudinal components respectively. The compari- 
son between the behavior of the transverse SVHF in the 
nematic (Fig. [5k.) and in the columnar (Fig. [5p, [5b, 
and [5b!) phase reveals a drastic change in the dynamics. 
The SVHF in the nematic phase is a monotonic function 
which broadens with time. By entering the columnar 
phase one observes the appearance of peaks which corre- 
spond to the positions of the hexagonal lattice in the xy 
plane. As expected, after a fixed time interval the num- 
ber and height of the peaks depend on the packing frac- 
tion of the system, so that by increasing the pressure the 
number of peaks decreases due to higher energetic bar- 
riers. These results confirm what was already observed 
for the smectic phase in experiments [34[ , simulation (42j 
and theory [39| , i.e. the partial translational symmetry 
breaking in a liquid crystal gives rise to a non-Gaussian 
quasi-quantized diffusion related to a hopping-type dy- 
namics. 

In Ref. [39l ] it was shown that, in order to accurately 
describe the dynamics of a liquid crystal system, it is not 
sufficient to take into account the permanent barriers due 
to the long-range structure, but also the transient caging 



FIG. 6. Longitudinal component of the self-part of the van 
Hove function Gl(z,t) as a function of z/D for t/r = 20 
and pressure P* = 2.5 (a), P* = 3.0 (b), P* = 3.5 (c) and 
P* — 4.0 (d). The curves refer to Gaussian fits over points 
near the origin (solid lines) and the tails (dashed lines). 



effect given by the surrounding particles. In this sense, 
the local fluid structure can affect the diffusion by deter- 
mining dynamical heterogeneities which make the system 
deviate from Gaussianity. A careful analysis on the lon- 
gitudinal component of the SVHF in Fig. [6] shows that 
this is indeed the case for the present columnar system. 
In fact, if along the z axis the diffusion was Gaussian, 
it would be possible to fit the points in Fig. [5] with a 
single Gaussian function. On the contrary, by perform- 
ing this fit on different intervals on the z axis, i.e. in 
the region near the origin (solid curve in figure) and the 
tails (dashed curve), one can observe that two different 
curves arc obtained. Although in the present system the 
deviations between the two curves are small, this behav- 
ior manifests interesting resemblances with the hetero- 
geneous dynamics of some amorphous systems, such as 
supercooled liquids and gels, where the two-Gaussian fit- 
ting is used to distinguish between "slow" and "fast" par- 
ticles [111 [59|]. In this sense, one can affirm that in the 
longitudinal direction the diffusion can be regarded as 
that of a dilute supercooled liquid, more than a normal 
liquid. The effect described so far should not be con- 
sidered as strictly due to the columnar structure of the 
system, since analogous deviations are also observed in 
Fig. UK for the nematic phase. Instead, the high packing 
fraction causes these small discrepancies from Gaussian 
diffusion. 

A description of the transient caging regime due to 
the nearest-neighbor (solvation) shell around each parti- 
cle can be given in terms of the distinct-part of the van 
Hove function (DVHF) denned in Eq. and © and 
reported in Fig. [7]for t/r = 0.02, 2 and 20. According to 
the definition given in Eq. ©, at time t = the DVHF 
coincides with the pair distribution function, and it is 




FIG. 7. (Color online) Transverse Gj[{R,t) (a-d) and longi- 
tudinal G\(z,t) (e-h) component of the distinct-part of the 
van Hove function evaluated for the same system and state- 
points as in Fig. [5] and [6] and at t/r — 0.02 (blue dotted 
line), t/r = 2 (green dashed line), t/r = 20 (red solid line) at 
P* = 2.5, 3.0, 3.5 and 4.0 (from top to bottom). 



thus characterized by a region around the origin where 
its value is equal to zero due to the excluded volume in- 
teraction. On the other hand, in the limit t — > 00 the 
DVHF is expected to be a constant in a translationally 
homogeneous system due to the decay of the positional 
correlations; this is not the case in presence of transla- 
tional order, since the mutual position of two particles at 
different times is influenced by the permanent long-range 
structure of the whole system. At t/r = 0.02 a region 
around the origin where the DVHF is close to zero sug- 
gests that each particle is still rattling around its initial 
position. Beyond this region a series of peaks indicate 
the preferential positions of the particles with respect to 
the one placed at the origin at the initial time. In the 
nematic phase at t/r = 0.02 one can recognize in both 
the transverse (Fig. [7^) and longitudinal (Fig. [7k) com- 
ponents the liquid-like structure of the system, where the 
lack of long-range order is testified by the rapid decay of 
the peaks by moving away from the origin. One should 
also notice the speed at which the gap region around the 
origin is filled, giving rise to a DVHF almost constant 
already at t/r = 2. Therefore, during this time interval 
a given particle i can escape the trapping cage formed by 
its nearest neighbors j, and the space originally occupied 
by i will be filled by one of the j particles. The situa- 
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FIG. 9. Non-Gaussian parameter {a,2,xy){t) of the diffusion in 
the plane perpendicular to the nematic director as a function 
of tjr for the same system and statepoints as in Fig. 08] In 
the inset we compare at P* = 3.0 the non-Gaussian parame- 
ter averaged over the concentrations (open circles) with those 
relative to species 1 (dashed line), to species 2 (dotted line) 
and of the whole system (solid line). 



FIG. 8. Mean square displacement in the direction perpendic- 
ular (AR 2 ) (a) and parallel (Az 2 ) (b) to the nematic director 
in units of D 2 as a function of t/r for the same system and 
statepoints as in Fig. 15171 



tion changes appreciably when we pass to the columnar 
phase, where the long-range modulations in the trans- 
verse component of the DVHF (Fig. [7p-d) indicate the 
presence of a permanent structure. On the other hand, 
the longitudinal component (Fig. [7f-h) does not display 
any dramatic change in shape, but from its time evolu- 
tion one can observe that the time a particle needs to 
leave its initial position increases considerably. In fact, 
whereas the relaxation times of the longitudinal DVHF 
in the three systems manifesting columnar order are com- 
parable (Fig. Ef-h), one can notice a faster relaxation in 
the nematic phase (Fig. [7J:), which cannot be due to the 
difference in packing fraction exclusively. This seems to 
suggest that the inhomogeneous structure and the result- 
ing dynamics in the transverse plane appreciably affect 
the dynamics in the longitudinal direction. We argue 
that the higher in-plane mobility of the nematic phase 
with respect to the columnar affects the mobility along 
the nematic director albeit only slightly. This is coherent 
with the results of Ref. [39| , where a coupling between 
transverse and longitudinal diffusion in the smcctic phase 
was pointed out. 

An alternative way to analyze dynamical hetero- 
geneities is to look for deviations from linearity of the 
mean square displacement. The effect of local cage trap- 
ping in systems close to dynamical arrest and the pres- 
ence of permanent long-range inhomogencitics as in liq- 
uid crystals manifest themselves in a region at interme- 



diate times where the dynamics is strongly sub-diffusive. 
In Fig. [5] we show the MSD both in the xy plane and 
in the z direction. In the plane perpendicular to the ne- 
matic director (Fig. [8^,) one can appreciate the almost 
linear trend of the MSD in the nematic phase, whereas 
by increasing the pressure and going to the columnar 
phase a plateau region appears, manifesting the devel- 
opment of a heterogeneous dynamics. These deviations 
from linearity are tightly related to the non-Gaussian be- 
havior of the self-part of the van Hove function, and can 
be quantitatively estimated by the non-Gaussian param- 
eter defined in Eq. (|9|). In Fig. [9] we report the NGP 
in the xy plane averaged over the species concentrations 
as described in Eq. (fTU)) . This parameter remains close 
to zero in the nematic phase, while it displays a peak 
at intermediate times in the columnar phase indicating 
deviations from Gaussianity. On the other hand, along 
the z direction (not shown here) the NGP does not devi- 
ate significantly from zero. The choice of calculating the 
NGP for the whole system by averaging over the value it 
assumes for the two species separately allows to take into 
account just the effects related to the long-range struc- 
ture of the system. For the sake of completeness we show 
in the inset of Fig. |H] a comparison between the NGP at 
P* = 3.0 for each species, their weighted average and 
that corresponding to the whole system. As expected, 
the operation of average does not affect significantly the 
position of the peak, but decreases only the peak height, 
suggesting that in this way the non- Gaussianity due to 
particle size difference is subtracted. One should notice 
that this particular treatment is not necessary for the rest 
of the physical properties measured in this paper, since 
with the exception of the NGP they result to be linear in 
the particle species, i.e. their value for the system as a 
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whole corresponds to a weighted average over the species. 

The dynamic inhomogeneities as captured by both the 
MSD and the NGP allow to identify three different time 
intervals. At short times the MSD follows the usual lin- 
ear trend and the NGP maintains a value close to zero, 
which means that particles diffuse freely since they do not 
feel yet the trapping cage due to the surrounding parti- 
cles. At intermediate times the MSD becomes strongly 
sub-diffusive and the NGP is characterized by a mono- 
tonic growth, thus meaning that the free diffusion is in- 
hibited by the columnar structure of the fluid. At this 
stage one can distinguish between particles which still 
rattle around the position of their column and others 
which succeeded in overcoming the energetic barrier and 
jumped to another column. The end of the subdiffusivc 
plateau and the return to the linear trend of the MSD 
roughly correspond to the peak of the NGP, which starts 
decreasing monotonically to zero indicating the end of 
the caging regime, i.e. most of the particles succeeded in 
leaving their initial column. A deeper inspection on the 
pressure dependence of the NGP shows that the degree 
of non-Gaussianity, i.e. the height of the peak, and the 
duration of the caging regime, i.e. the position of the 
peak, increase with packing fraction. This fact can be 
explained by considering that the cage escape is related 
to a rearrangement of the surrounding particles, which 
becomes slower at higher packing fraction as it involves 
more of them. Furthermore, the small deviations from 
linearity in the MSD in the direction parallel to the ne- 
matic director confirm the presence of a weakly hetero- 
geneous dynamics, as already pointed out by analyzing 
the self-part of the van Hove function in this direction. 

Finally, the structural relaxation of the system is an- 
alyzed in terms of the self-part of the intermediate scat- 
tering function (SISF) defined in Eq. (TTTj) . Whereas 
along the z direction the relaxation is characterized by a 
single step decay at each pressure (Fig. [TUb). a plateau 
region, which characterizes the relaxation in the xy plane 
at intermediate times, develops in the columnar phase. 
This plateau, whose value increases with pressure, indi- 
cates the time extension of the cage regime and is ex- 
pected to divide a short-time decay (/3-relaxation) from 
a long-time one (a-relaxation). As previously observed 
in recent work on smectic liquid crystals |4lM43j and in 
out-of-equilibrium supercooled liquids j47j, the SISF de- 
cays likely to zero at long times, indicating the loss of 
density auto-correlations. This kind of behavior was de- 
scribed for the smectic phase in Refs. (4l|-|43|. where 
the a-relaxation decay was fitted by a stretched expo- 
nential function of the form ex-p[(t / t r )°] with /3 ~ 0.6 
and t r the characteristic relaxation time. In the present 
simulations we did not observe any a-relaxation as the 
relaxation time of the systems exceeds probably our simu- 
lation time. On the other hand, from the data available a 
close accordance with the features of the structural relax- 
ation of the smectic phase can be observed. In particular, 
the /3-relaxation in the xy plane is reasonably described 
by an exponential decay, as expected for simple liquids, 
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FIG. 10. Self-part of the intermediate scattering function 
F s ,xy{t) and F SiZ (t) evaluated at wave vectors corresponding 
to the first peak of the structure factor in the plane perpen- 
dicular to the nematic director (a) and along it (b). The data 
correspond to the same system and statepoints as in Fig. 15191 
The solid lines are fits (see text). 



due to the lack of interactions of the particles with the 
nearest neighbors at small times. Also the relaxation 
along the z axis resembles accurately the relaxation of 
the smectic phase inside the smectic layers. In fact, in 
both these cases the SISF depends weakly on the pres- 
sure and it is characterized by an exponential decay at 
small times, which eventually becomes a stretched ex- 
ponential with /? ~ 0.6. In this sense we can confirm 
what the authors observed in Ref. [4l|, i.e. the relax- 
ation of a liquid crystal in the direction(s) in which the 
system is homogeneous is closer to that of a low-density 
supercooled liquid, than a simple liquid, where instead 
an exponential relaxation is to be expected. 



IV. CONCLUSIONS 

In summary, we used Monte Carlo simulations to ana- 
lyze for the first time the presence of dynamical hetero- 
geneities in a columnar liquid crystal of perfectly aligned 
hard spherocylindcrs. The long-range hexagonal order 
in the plane perpendicular to the nematic director de- 
termines an effective mean-field potential, whose effect 
is to maintain particles inside a column preventing them 
to occupy a position in between the columns. In analogy 
with previous analyses on the smectic phase, the height of 
the energetic barriers of this effective potential increases 
with the packing fraction and the particle anisotropy. As 
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a consequence, in the xy plane the dynamics of a rod 
is characterized by a quasi-quantized behavior in which 
particles rattle around the position of the column and 
jump to another column only when the configuration of 
the surrounding particles allows it. 

The rattling-and-jumping dynamics in the in-plane 
evolution of the system gives rise to three different time 
regimes. At very short times, the particles diffuse al- 
most freely because they do not feel yet the presence 
of the trapping cage formed by their surrounding near- 
est neighbors. At this stage the behavior of the system is 
that typical of a simple fluid, characterized by a Gaussian 
distribution of displacements, a linear mean square dis- 
placement and a fast exponential structural relaxation. 
A second stage starts when particles begin experiencing 
the cage due to the long-range structure of the system, in 
such a way that the diffusion results to be inhibited and 
only occasionally a column-to-column jump takes place 
and is made possible by the instantaneous configuration 
of the system. As a result, the mean square displace- 
ment as well as the self-intermediate scattering function 
develop a plateau, which testifies the slowing down of 
the dynamics and whose time extension increases with 
packing fraction. On the other hand, the distribution of 
displacements show marked deviations from Gaussianity 
due to the appearance of peaks which correspond to the 
lattice positions in the plane. Nonetheless, after longer 
time intervals the number of "fast" particles, which suc- 
ceeded in overcoming the energetic barrier, increases with 
respect to the "slow" ones. Consequently, when most of 
the particles succeeded in leaving their initial column, a 
second diffusive regime starts, indicating the end of the 
cage regime. 

We observed interesting analogies with the dynamics 



in smectic phases by considering the in-column dynamics. 
In fact, along the direction defined by the nematic direc- 
tor, the system does not develop any long-range order, 
and it is thus expected to behave like a liquid. On the 
other hand, we noticed interesting, although slight, de- 
viations from Gaussian diffusion both in the distribution 
of displacements and in the mean square displacement. 
As far as the structural relaxation is concerned, this fact 
is testified by a self-intermediate scattering function well 
approximated by a stretched-exponential, as it happens 
in dense liquids. In this sense we confirm previous studies 
on the smectic phase, that along the direction in which a 
liquid crystal does not develop any long-range order the 
dynamics is similar to a dense liquid. These results are 
to be compared with recent experiments on the colum- 
nar phase of a suspension of fd virus particles |44| , where 
huge discrepancies from Gaussianity were observed along 
the nematic director. We argue that the higher length- 
to-diameter ratio, the flexibility or the charge of the rods 
could account for a more pronounced non-Gaussian dif- 
fusive behavior than what was observed in the present 
study. 
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